home *** CD-ROM | disk | FTP | other *** search
/ PC Media 23 / PC MEDIA CD23.iso / share / prog / newmat / newmatnl.cpp < prev    next >
C/C++ Source or Header  |  1995-01-11  |  6KB  |  184 lines
  1. //$$ newmatnl.cpp         Non-linear optimisation
  2.  
  3. // Copyright (C) 1993,4,5: R B Davies
  4.  
  5.  
  6. #define WANT_MATH
  7. #define WANT_STREAM
  8.  
  9. #include "newmatap.h"
  10. #include "newmatnl.h"
  11.  
  12.  
  13.  
  14.  
  15. void FindMaximum2::Fit(ColumnVector& Theta, int n_it)
  16. {
  17.    Tracer tr("FindMaximum2::Fit");
  18.    Real z,w,x,x2,g,l1,l2,l3,d1,d2,d3;
  19.    ColumnVector Theta1, Theta2, Theta3;
  20.    int np = Theta.Nrows();
  21.    ColumnVector H1(np), H3, HP(np), K, K1(np);
  22.    Boolean oorg, conv;
  23.    int counter = 0;
  24.  
  25.    {                              // subblock so won't have labels in same
  26.                                   // block as destructors.
  27.  
  28.                                   // I know it is supposed to be evil to
  29.                                   // use "goto". Nevertheless this seems
  30.                                   // to be the best way of coding the
  31.                                   // algorithm, even in C++.
  32.  
  33.       Theta1 = Theta; HP = 0.0; g = 0.0;
  34.  
  35.    Start:
  36.       tr.ReName("FindMaximum2::Fit/Start");
  37.       Value(Theta1, TRUE, l1, oorg); if (oorg) Throw(DataException(Theta));
  38.  
  39.    Restart:
  40.       tr.ReName("FindMaximum2::Fit/ReStart");
  41.       conv = NextPoint(H1, d1); if (conv) goto Convergence;
  42.       if (counter++ > n_it) goto Fail;
  43.  
  44.       z = 1.0 / sqrt(d1);
  45.       H3 = H1 * z; K = (H3 - HP) * g; HP = H3;
  46.       g = 0.0;                     // de-activate to use curved projection
  47.       if (g==0.0) K1 = 0.0; else K1 = K * .2 + K1 * .6;
  48.       // (K - K1) * alpha + K1 * (1 - alpha) = K * alpha + K1 * (1 - 2 * alpha)
  49.       K = K1 * d1; g = z;
  50.  
  51.    Continue:
  52.       tr.ReName("FindMaximum2::Fit/Continue");
  53.       Theta2 = Theta1 + H1 + K;
  54.       Value(Theta2, FALSE, l2, oorg);
  55.       if (counter++ > n_it) goto Fail;
  56.       if (oorg)
  57.          { H1 = H1 * 0.5; K = K * 0.25; d1 *= 0.5; g *= 2.0; goto Continue; }
  58.       d2 = LastDerivative(H1 + K * 2.0);
  59.  
  60.    Interpolate:
  61.       tr.ReName("FindMaximum2::Fit/Interpolate");
  62.       z = d1 + d2 - 3.0 * (l2 - l1);
  63.       w = z * z - d1 * d2;
  64.       if (w < 0.0) goto Extrapolate;
  65.       w = z + sqrt(w);
  66.       if (1.5 * w + d1 < 0.0) goto Extrapolate;
  67.       if (d2 > 0.0 && l2 > l1 && w > 0.0) goto Extrapolate;
  68.       x = d1 / (w + d1); x2 = x * x; g /= x;
  69.       Theta3 = Theta1 + H1 * x + K * x2;
  70.       Value(Theta3, TRUE, l3, oorg);
  71.       if (counter++ > n_it) goto Fail;
  72.       if (oorg)
  73.       {
  74.          if (x <= 1.0)
  75.             { x *= 0.5; x2 = x*x; g *= 2.0; d1 *= x; H1 = H1 * x; K = K * x2; }
  76.          else
  77.          {
  78.             x = 0.5 * (x-1.0); x2 = x*x; Theta1 = Theta2;
  79.             H1 = (H1 + K * 2.0) * x;
  80.             K = K * x2; g = 0.0; d1 = x * d2; l1 = l2;
  81.          }
  82.          goto Continue;
  83.       }
  84.  
  85.       if (l3 >= l1 && l3 >= l2) { Theta1 = Theta3; l1 = l3; goto Restart; }
  86.  
  87.       d3 = LastDerivative(H1 + K * 2.0);
  88.       if (l1 > l2)
  89.          { H1 = H1 * x; K = K * x2; Theta2 = Theta3; d1 = d1*x; d2 = d3*x; }
  90.       else
  91.       {
  92.          Theta1 = Theta2; Theta2 = Theta3;
  93.          x -= 1.0; x2 = x*x; g = 0.0; H1 = (H1 + K * 2.0) * x;
  94.          K = K * x2; l1 = l2; l2 = l3; d1 = x*d2; d2 = x*d3;
  95.          if (d1 <= 0.0) goto Start;
  96.       }
  97.       goto Interpolate;
  98.  
  99.    Extrapolate:
  100.       tr.ReName("FindMaximum2::Fit/Extrapolate");
  101.       Theta1 = Theta2; g = 0.0; K = K * 4.0; H1 = (H1 * 2.0 + K);
  102.       d1 = 2.0 * d2; l1 = l2;
  103.       goto Continue;
  104.  
  105.    Fail:
  106.       Throw(ConvergenceException(Theta));
  107.  
  108.    Convergence:
  109.       Theta = Theta1;
  110.  
  111.    }
  112. }
  113.  
  114.  
  115.  
  116. void NonLinearLeastSquares::Value
  117.    (const ColumnVector& Parameters, Boolean, Real& v, Boolean& oorg)
  118. {
  119.    Tracer tr("NonLinearLeastSquares::Value");
  120.    Y.ReDimension(n_obs); X.ReDimension(n_obs,n_param);
  121.    // put the fitted values in Y, the derivatives in X.
  122.    Pred.Set(Parameters);
  123.    if (!Pred.IsValid()) { oorg=TRUE; return; }
  124.    for (int i=1; i<=n_obs; i++)
  125.    {
  126.       Y(i) = Pred(i);
  127.       X.Row(i) = Pred.Derivatives();
  128.    }
  129.    if (!Pred.IsValid()) { oorg=TRUE; return; }  // check afterwards as well
  130.    Y = *DataPointer - Y; Real ssq = Y.SumSquare();
  131.    errorvar =  ssq / (n_obs - n_param);
  132.    cout << "\n" << setw(15) << setprecision(10) << " " << errorvar;
  133.    Derivs = Y.t() * X;          // get the derivative and stash it
  134.    oorg = FALSE; v = -0.5 * ssq;
  135. }
  136.  
  137. Boolean NonLinearLeastSquares::NextPoint(ColumnVector& Adj, Real& test)
  138. {
  139.    Tracer tr("NonLinearLeastSquares::NextPoint");
  140.    QRZ(X, U); QRZ(X, Y, M);     // do the QR decomposition
  141.    test = M.SumSquare();
  142.    cout << " " << setw(15) << setprecision(10)
  143.       << test << " " << Y.SumSquare() / (n_obs - n_param);
  144.    Adj = U.i() * M;
  145.    if (test < errorvar * criterion) return TRUE;
  146.    else return FALSE;
  147. }
  148.  
  149. Real NonLinearLeastSquares::LastDerivative(const ColumnVector& H)
  150. { return (Derivs * H).AsScalar(); }
  151.  
  152. void NonLinearLeastSquares::Fit(const ColumnVector& Data,
  153.    ColumnVector& Parameters)
  154. {
  155.    Tracer tr("NonLinearLeastSquares::Fit");
  156.    n_param = Parameters.Nrows(); n_obs = Data.Nrows();
  157.    DataPointer = &Data;
  158.    FindMaximum2::Fit(Parameters, Lim);
  159.    cout << "\nConverged\n";
  160. }
  161.  
  162. void NonLinearLeastSquares::MakeCovariance()
  163. {
  164.    if (Covariance.Nrows()==0)
  165.    {
  166.       UpperTriangularMatrix UI = U.i();
  167.       Covariance << UI * UI.t() * errorvar;
  168.       SE << Covariance;                 // get diagonals
  169.       for (int i = 1; i<=n_param; i++) SE(i) = sqrt(SE(i));
  170.    }
  171. }
  172.  
  173. void NonLinearLeastSquares::GetStandardErrors(ColumnVector& SEX)
  174.    { MakeCovariance(); SEX = SE.AsColumn(); }
  175.  
  176. void NonLinearLeastSquares::GetCorrelations(SymmetricMatrix& Corr)
  177.    { MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); }
  178.  
  179. void NonLinearLeastSquares::GetHatDiagonal(DiagonalMatrix& Hat) const
  180. {
  181.    Hat.ReDimension(n_obs);
  182.    for (int i = 1; i<=n_obs; i++) Hat(i) = X.Row(i).SumSquare();
  183. }
  184.